Environmental vs. demographic variability in two-species predator-prey models 
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We investigate the competing effects and relative importance of intrinsic demographic and envi- 
ronmental variability on the evolutionary dynamics of a stochastic two-species Lotka-Volterra model 
by means of Monte Carlo simulations on a two-dimensional lattice. Individuals are assigned inher- 
itable predation efficiencies; quenched randomness in the spatially varying reaction rates serves as 
environmental noise. We find that environmental variability enhances the population densities of 
both predators and prey while demographic variability leads to essentially neutral optimization. 

PACS numbers: 87.23.Cc, 05.40.-a, 87.18.Tt 



The mathematical modeling of species interactions 
continues to be a central issue in population ecology [H-13| ■ 
Several simple models have been proposed, investigated, 
and sometimes realized under laboratory conditions. Yet 
more realistic and thus biologically more relevant model 
variants obviously have to include both external spatial 
disorder in the reaction rates to account for varying envi- 
ronmental conditions and intrinsic demographic hetero- 
geneity stemming from trait variability in individuals. 
While we addressed the former in a recent study [5[ , our 
goal in this letter is to investigate the interplay between 
quenched spatial rate disorder and additional variability 
of individuals' reaction rates, as well as intriguing evolu- 
tionary co-optimization within interacting populations. 

We focus on the Lotka-Volterra (LV) predator-prey 
model owing to its simplicity and because its basic fea- 
tures are well-understood. It was first introduced to 
study fish populations in the Adriatic sea and chemi- 
cal oscillations 0, 0]. While the original deterministic 
LV (mean-field) equations yield neutral cycles and hence 
persistent nonlinear oscillations around a marginal fixed 
point 1], in stochastic implementations this species co- 
existence fixed point becomes stable and is approached 
very slowly through damped oscillations Spatially 
extended stochastic versions of the LV model yield strik- 
ing dynamical patterns and emergent inter-species cor- 
relations 15l421| which may be utilized to quantitatively 
assess the response to external or internal changes. Pop- 
ulation stability can be measured via the extinction time 
in small systems, where the stochastic kinetics ultimately 
reaches an absorbing zero-particle state 
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In our study of the effects of environmental rate vari- 
ability in the LV model, we found a remarkable increase 
of the asymptotic population densities of both species 
with enhanced quenched spatial disorder, i.e., predation 
rates that are fixed to different lattice sites [5]. Yet 
the observed erratic population oscillations and relax- 
ation towards the (quasi-)steady state occur on the time 
scale of many generations; for real biological systems, 
one therefore needs to address Darwinian evolutionary 
adaptation of individuals' traits. Consequently, we in- 
troduce fundamentally novel features by endowing in- 



dividual predator and prey particles with randomly se- 
lected rates, and investigate whether and how optimiza- 
tion within each species due to imperfect efficiency inher- 
itance (mimicking random mutations) further reinforces 
the total population's stability and fitness. Dynamical 
coevolution of interacting species is a crucial feature of 
adapting ecological systems and has been studied experi- 
mentally 23, l24[ as well as theoretically 25-28]. Combin- 
ing quenched spatial with individual, evolving rate dis- 
tributions allows us to quantitatively assess the relative 
importance of environmental vs. demographic, inheri- 
table variabilities in a nonlinear competing two-species 
predator-prey system. 

We find that both environmental and individual-based 
variabilities combined with random mutations produce 
a marked enhancement of the quasi-stationary densities 
of both species, thus considerably extending our earlier 
conclusions for purely environmental randomness In 
addition, individual variability stabilizes both predator 
and populations against extinction. Remarkably, the op- 
timization of predation and evasion capabilities of either 
species turns out to be essentially neutral in the popu- 
lation densities; in contrast to genetic drift models [29j, 
our nonlinear model does not lead to trait fixation. 

We consider a spatially extended version of the LV 
model consisting of two particles species. The "predator" 
species is subject to spontaneous decay A — > with rate 
/i, while the "prey" B reproduce (asexually): B — > 2B 
with rate a. Different particles interact on-site with a 
non-uniform predation rate A, whereupon a prey is re- 
moved and replaced with a predator: A + B — > 2A. The 
prey birth and predator death rates both remain fixed 
at a uniform value <j = fJ, = 0.5 for all particles and 
lattice sites, whereas the predation rates are allowed to 
vary between different positions and participating par- 
ticles (see below). Particles exist on a two-dimensional 
square lattice with 128 x 128 sites and periodic bound- 
ary conditions. (We could not find significant finite-size 
effects already at this lattice size.) Both species perform 
unbiased random walks via nearest-neighbor hopping oc- 
curring with probability one, hence all rates are to be 
understood relative to the diffusivity D. Reactions occur 
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on-site, assuming infinite local carrying capacities, imply- 
ing that the growth of the population on any single site is 
essentially unrestricted (with a safety limit of rij < 1000 
per lattice site i that is never reached with the param- 
eters investigated in the present study). The predator 
extinction transition occurring in model variants with re- 
stricted site occupation is thus absent here [12, 2(J 21 1. 
The initial population distribution of both predator and 
prey particles is chosen randomly with a mean density 
PA,i = Pba = 1- The simulation proceeds via random se- 
quential updates, with one Monte Carlo step being com- 
pleted when on average each particle in the simulation 
has moved and had a chance to react [30j ] . 

In order to model variability of individuals and trait in- 
heritance, each particle carries a predation efficacy prop- 
erty T] £ [0, 1], determined during the particle's creation 
and providing a coarse-grained characterization of the 
combined efficacies of its genetic heritage (genes) and its 
learned strategies (memes). An offspring's position in 
efficiency space will thus be near its parent's location 
but subject to random changes (mutations in the case of 
genes, adaptations of strategies in the case of memes), 
thereby suggesting the use of a normalized Gaussian dis- 
tribution centered at the parent's efficiency value r\p 
(truncated to the interval [0, 1] accessible to a reaction 
probability) to assign an efficiency value 770 to the off- 
spring during reproduction. The standard deviation wp 
of the Gaussian function constitutes a model parame- 
ter and corresponds to the average severity of mutations 
from one generation to the next. Note that the efficiency 
assigned to a particle 77 is different from the traditional 
genetic fitness, which is defined as the average number of 
offspring produced by a genome. It represents a meso- 
scopic continuous stochastic variable, as opposed to a 
genetic description employing naturally discrete values. 

Since we wish to address the distinctions between in- 
ternal and spatial randomness, we introduce in addition 
environmental variability by assigning a spatial preda- 
tion efficacy value 775 to each lattice site, drawn from 
a normalized Gaussian distribution with fixed mean 0.5 
and standard deviation ws, truncated to [0,1], and set 
to be fixed in time [3lJ. The ensuing predation rate A is 
a random variable as well, namely a function of both the 
spatial efficiency at the lattice site the reaction occurs on 
and the two individual predation efficacies of the partic- 
ipating predator and prey particles. We finally define a 
model parameter £ that describes the relative importance 
of the spatial over individual efficacies: 



A = Ct?s + (1-C) (va + Vb)/2. 



(1) 



Over many generations, species optimize their pre- 
dation / evasion efficiency through evolving their ef- 
ficacy distributions by means of random inheritance. 
Figure HJa) shows the population density histograms 
Pa/b{v) f° r a representative case with moderate inher- 
itance variability wp = 0.1. The initial population of 




FIG. 1: (Color online.) (a) The predator (red) and prey (blue) 
population densities in the (quasi-)steady state as functions 
of the predation efficiency, averaged over 10 4 Monte Carlo 
simulation runs, optimize towards high rj A w 0.71 and low 
rj B 0.26 mean values, respectively. The standard deviation 
of the Gaussian inheritance distribution is wp = 0.1 and the 
spatial influence factor £ = 0. In contrast to genetic drift 
models, both densities do not fixate at extreme predation ef- 
ficiencies (77 = 0, 1). The dashed lines show the theoretical 
prediction from an effective stochastic mean-field model, (b) 
For a flat inheritance distribution (wp — 00), the predators 
experience no selection bias while the prey population is pref- 
erentially selected towards low predation efficacies. 



predator and prey particles had an assigned predation ef- 
ficacy of tia/b = 0.5. We have carefully checked that the 
final (quasi-)steady state population distribution does 
not depend on the (uncorrelated) initial conditions (ex- 
cept for those rare simulation runs when either the prey 
or predator population went extinct). The predator pop- 
ulation maximum moves towards higher mean efficacy 
whereas the prey population, for which lower values of 
the predation efficiency are preferable, tends towards a 
lower average. Predators with a slightly higher efficiency 
value are more successful at predation and thereby re- 
produce more often. Hence their improved predation ca- 
pability is passed on to subsequent generations with a 
higher frequency, driving the overall predator population 
toward higher efficiency values. Similarly, prey particles 
with a lower predation efficiency are better at evasion 
and thus survive longer. This gives them the chance to 
reproduce at a higher rate, driving the prey population 
towards low mean efficacy. In the extreme situation of 
completely random assignment of predation efficiencies, 
where no correlations between the corresponding values 
for parents and offspring are implemented (equivalent to 
a uniform inheritance distribution with wp — 00), we al- 
ready see a strong tendency towards low efficacies for the 
prey species; see Fig. [T](b). This feature is explained by 
the bias in the predation rule that favors selection of prey 
particles with higher efficiency values. For predators no 
such bias exists; hence their population distribution in 
efficiency space remains flat. Spatial fluctuations modify 
the results quantitatively, but not qualitatively, whence 
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FIG. 2: (Color online.) Mean extinction time t e of a small 
system of 10 x 10 lattice sites as a function of individual vari- 
ability wp and for ( = 0. Extinction here is defined as the 
event when either species goes extinct. With increasing vari- 
ability the system on average takes longer to go extinct. 

we observe a slightly more pronounced effect in the prey 
distribution in our studies of non-spatial systems. 

In order to analytically verify our data we write 
down the mean-field equations for a well-mixed zero- 
dimensional LV system with individual variability. We 
consider the number of particles of species A and B as a 
function of predation efficiency. The predation efficiency 
77 is a bounded quantity in the interval [0, 1] assigned 
to particles at their time of creation. We discretize this 
interval into N bins with a spacing of A77 = 1/N and 
midpoint rji = (i + 1/2) /N, and denote the number of 
A and B particles in a bin i respectively as a% and bi 
(i = 0, . . . , N — 1). To model individual variability we in- 
troduce the probability = f(T)i,r)j) for a parent with 
efficiency rji to produce offspring with efficiency rjj . The 
predation rate is a function of the efficacies of the preda- 
tor A and prey B participating in the predation reaction: 
Kj = (Vi + 7 lj)/2 [essentially the discretized Eq. ((T|) with 
C = 0]. Thus we arrive at the coupled mean- field rate 
equations for the case of purely individual variability: 

hi = -fiat + J2-j Hk ^jkfikakbj , (2) 

b i = a Hk fikh ~ J2j ^ij a jbi ■ (3) 

The steady-state densities are obtained by setting the 
time derivatives to zero, yielding expressions for <Zj and 
bi that can be solved iteratively for any inheritance prob- 
ability distribution /, as shown (dashed) in Fig.QJa). 

In the special case of a uniform probability distribu- 
tion (implying the absence of any correlation between 
the predation efficiencies of parent and offspring parti- 
cles) fij = 1/N, the steady-state densities, defined by 
PA.i = ai/J2j a i and PB.i = h/ J2j bj, can be obtained 
exactly. The predator population acquires a constant (77- 
independent) value of pa = 1/N, whereas the prey popu- 
lation decreases with increasing 77 as pb(ti) — N f n3 \+2-q ■ 
This result is exactly mirrored by our zero-dimensional 
Monte Carlo simulations. In spatially extended systems 
fluctuations modify the density distributions, leading to 
a prey density dependence that is slightly less steep as a 
function of the predation efficiency, see Fig. [Tfb) . Corre- 
lation effects not captured by mean-field theory are evi- 
dently strongest at the distribution maxima. 




FIG. 3: (Color online.) (a) The (quasi-)steady state predator 
density pA as a function of the spatial and individual variabil- 
ity ws and wp for £ = 0.3. The black line represents the slice 
of equal spatial and individual variability ws = Wp from the 
minimum in (b). (b) The predator density shows a consis- 
tent increase for all values of the spatial variability influence 
£ as function of equal variabilities w = ws — wp over a sys- 
tem with zero variability. A remarkable minimum is observed 
near £ = 0.3 (black line), (c) The standard deviation of the 
predation rate o\, calculated via error propagation from the 
spatial and individual predation efficiency distributions. 



We collected extinction time histograms for small sys- 
tems (lattice size 10 x 10 sites) to determine the influ- 
ence of individual variability on the stability of the pop- 
ulation. In finite stochastic systems with an absorbing 
state (here, predator extinction), fluctuations will even- 
tually drive the system into the absorbing state. Figure [2] 
demonstrates that the mean extinction time is enhanced 
by a factor up to w 4.5 by individual variability, render- 
ing the system markedly more robust against extinction. 

To quantify the influence of variability, and in particu- 
lar the distinction between individual (internal) variabil- 
ity and spatial environmental randomness, we measured 
its impact on the (quasi-)steady state particle density for 
both species. Figure [3^a) displays the relative change of 
the predator density pa{ws, wp) over the zero- variability 
case as a function of ws and wp for £ = 0.3. Both types 
of variability contribute additively and positively to the 
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density enhancement. Figure^b) shows the relative den- 
sity change as a function of w = ws — wp and £ [32j |. 
The prey density shows the same quantitative behavior 
for all parameter ranges. Hence we observe a significant 
increase of the population densities of both species for 
higher variability not only for purely spatial = 1) ran- 
domness but also for individual variability (£ = 0). 
In contrast, the effect of spatial randomness in either the 
prey birth rate a or the predator death rate /x on the 
species densities stayed below a rather low value of 2%. 

The striking minimum in the density increase occur- 
ring near a spatial influence factor £ = 0.3 arises from the 
combined variabilities through the quenched randomness 
of the lattice and the emergent variability of the indi- 
vidual particles. We argue that the density increase is 
primarily a monotonic function of the variability in the 
predation rate A. Using the dependence of the preda- 
tion rate A on the spatial predation efficacy value rjs and 
the predation efficiencies of the participating particles t\a 
and t]b given in the text, the standard deviation of A is 
c \ = y/C 2 a' 2 s + (1 - C) 2 {o-'a + <J 2 B )/2. Due to the trun- 
cation of predation efficiency values to the range [0, 1], 
the effective standard deviation of the spatial predation 
efficacy is different from the environmental variability pa- 
rameter. Similarly, the standard deviation of the preda- 
tion efficiencies of individual particles have to be taken 
from simulation data. Figure 02c) shows the resulting 
standard deviation of A as a function of w and £ which 
is a measure of the effective combined variability. It re- 
flects the minimum in the density increase at £ rs 0.3. 
The data also emphasize that environmental variability 
has a more pronounced effect on the species densities as 
compared to demographic variability, since the density 
increase is disproportionally higher for £ — > 1. 

Surprisingly, we observe that low individual variability 
with weak or no spatial influence, i.e. < wp <C 1 and 
C = 0, yields the strongest species optimization with the 
maxima of the predator and prey populations closest to 
r] = 1 and ?y = respectively; see Fig. [Ha). But the en- 
hancement of the overall species densities in this regime is 
minute and tends to zero for small w p ; see the lower right 
corner of Fig. |3jb). The respective benefits of the up- / 
downward optimization of the predator / populations in 
terms of predation efficiency clearly almost cancel each 
other. Hence we conclude that predation efficiency opti- 
mization is essentially neutral and carries no benefit for 
either species in terms of their net population densities 
(at least in the context of our model), despite its vital ne- 
cessity to ensure the survival of coevolving species. This 
also reinforces our argument that the density enhance- 
ment is a function of rate variability only. 

The predator-predator and prey-prey correlation 
lengths Iaa and Ibb and typical predator-prey distance 
Iab, measured by extracting the (quasi-)steady state ex- 
ponential decay length and the position of the maxi- 
mum (in the case of Iab) from the correlation functions 




FIG. 4: (Color online.) (a) The correlation length Iaa from 
the predator-predator autocorrelation function. The prey- 
prey correlation length Ibb would essentially display the same 
shape as Iaa scaled by ~ 0.9. (b) The predator density re- 
laxation time T re iax toward the quasi-stationary state. 

C ap >(x) = (n ai+x n/3i) ~ p a pp with a,/3 — A,B (33||, de- 
crease for increasing variability w; see Fig. [4f a) . For 
C = 1 we reproduce the data from Ref. [5(, where we 
argued that the decrease in l a p indicated a more tightly 
clustered population around lattice sites with small spa- 
tial predation efficiency 775, leading to the observed en- 
hanced densities and higher amplitudes in the initial pop- 
ulation oscillations. Surprisingly we also see a (less pro- 
nounced) decrease of l a p, for £ — >• 0, indicating the ex- 
istence of spontaneously formed tight activity patches 
around clusters of highly optimized prey particles. To 
investigate the effect of the combined variability on the 
relaxation of the population densities we collected data 
on the characteristic decay time of the initial predator 
population oscillations by least-square fitting of an ex- 
ponentially decaying sinosiodal function to the predator 
species density time series; see Fig.[4|b). As expected, in- 
creasing disorder w induces a roughly threefold decrease 
in the purely spatial case (£ = 1) and about a twofold 
decrease in the individual variability case (£ = 0). 

In conclusion we performed an extensive numerical 
Monte Carlo simulation study to assess how external en- 
vironmental randomness and individual variability mod- 
ified through mutations during inheritance compete and 
affect the coevolutionary population dynamics of two co- 
existing species in a spatial stochastic LV model. The 
overall predator and prey densities are both enhanced by 
environmental variations, while evolutionary optimiza- 
tion within each species has an essentially neutral net 
effect. To better understand this evolutionary trait op- 
timization, we derived a mean-field model that qualita- 
tively reproduces our simulation results. In addition we 
find that increased individual variability stabilizes both 
populations against extinction. There are certainly other 
intriguing aspects pertaining to variability in ecologi- 
cal models that deserve further investigation, promising 
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amazingly rich features and crucial quantitative insight. 
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